/**************************************************************\
 **     Created originally by Jonathan Naylor, G4KLX.        **
 **     Later embellished by John Magliacane, KD2BD to       **
 **     detect and handle voids found in the SRTM data,      **
 **     SRTM-3 data in .BIL and .HGT format, and high        **
 **     resolution SRTM-1 one arc-second topography data.    **
 **************************************************************
 **                    Compile like this:                    **
 **      cc -Wall -O3 -s -lbz2 srtm2sdf.c -o srtm2sdf        **
 **              Last modification: 12-Mar-2009              **
\**************************************************************/ 

/*\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
Ported to Windows by VE3NCQ, Austin Wright, with several solutions from
a previous port by John McMellen. This program requires the proper Shuttle
Terrain data files for proper operation. See SPLAT documentation. This
program must be compiled to two different executables, one of which must be
srtm2sdf-hd.exe for use with 1-arc second terrain data. Put them in the same
folder as the downloaded and extracted .hgt files. A Windows GUI to batch
convert these files will be/is available from my website at www.ve3ncq.ca.
Consider copying the HD files in their own folder.
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <fcntl.h>
#include <unistd.h>
//#include <bzlib.h> removed BZ2 support by VE3NCQ to solve compile problems

#define BZBUFFER 65536

char	sdf_filename[30], sdf_path[255], replacement_flag, opened=0,
	hgt=0, bil=0;

int	srtm[3601][3601], usgs[1201][1201], max_north, max_west, n,
	min_north, min_west, merge=0, min_elevation, bzerror, ippd, mpi;

int ReadSRTM(char *filename)
{
	int x, y, infile, byte=0, bytes_read;
	unsigned char error, buffer[2];
	char north[3], west[4], *base=NULL, blw_filename[255];
	double cell_size, deg_north, deg_west;
	FILE *fd=NULL;

	if (strstr(filename, ".zip")!=NULL)
	{
		fprintf(stderr, "*** Error: \"%s\" must be uncompressed\n",filename);
		return -1;

	}

	if (strstr(filename, ".tgz")!=NULL)
	{
		fprintf(stderr, "*** Error: \"%s\" must be uncompressed\n",filename);
		return -1;

	}

	if ((strstr(filename, ".hgt")==NULL) && (strstr(filename, ".bil")==NULL))
	{
		fprintf(stderr, "*** Error: \"%s\" does not have the correct extension (.hgt or .bil)\n",filename);
		return -1;
	}

	if (strstr(filename, ".hgt")!=NULL)
		hgt=1;

	if (strstr(filename, ".bil")!=NULL)
		bil=1;

	base=strrchr(filename, '\\'); //changed to escaped "\" for windows by VE3NCQ

	if (base==NULL)
		base=filename;
	else
		base+=1;

	if (hgt)
	{
		/* We obtain coordinates from the base of the .HGT filename */

		north[0]=base[1];
		north[1]=base[2];
		north[2]=0;

		west[0]=base[4];
		west[1]=base[5];
		west[2]=base[6];
		west[3]=0;

		if ((base[0]!='N' && base[0]!='S') || (base[3]!='W' && base[3]!='E'))
		{
			fprintf(stderr, "*** Error: \"%s\" doesn't look like a valid .hgt SRTM filename.\n", filename);
			return -1;
		}

		max_west=atoi(west);

		if (base[3]=='E')
			max_west=360-max_west;

		min_west=max_west-1;

		if (max_west==360)
			max_west=0;

		if (base[0]=='N')
			min_north=atoi(north);
		else
			min_north=-atoi(north);

		max_north=min_north+1;
	}

	if (bil)
	{
		/* We obtain .BIL file coordinates
		   from the corresponding .BLW file */

		strncpy(blw_filename,filename,250);
		x=strlen(filename);

		if (x>3)
		{
			blw_filename[x-2]='l';
			blw_filename[x-1]='w';
			blw_filename[x]=0;

			fd=fopen(blw_filename,"rb");

			if (fd!=NULL)
			{
				n=fscanf(fd,"%lf",&cell_size);

				if ((cell_size<0.0008) || (cell_size>0.0009))
				{
					printf("\n*** .BIL file's cell size is incompatible with SPLAT!!\n");
					exit(1);
				}

				n=fscanf(fd,"%lf",&deg_west);
				n=fscanf(fd,"%lf",&deg_west);
				n=fscanf(fd,"%lf",&deg_west);

				n=fscanf(fd,"%lf",&deg_west);

				n=fscanf(fd,"%lf",&deg_north);

				fclose(fd);
			}

			min_north=(int)(deg_north);
			max_north=min_north+1;

			if (deg_west<0.0)
				deg_west=-deg_west;
			else
				deg_west=360.0-deg_west;

			min_west=(int)(deg_west);

			if (min_west==360)
				min_west=0;

			max_west=min_west+1;
		}
	}
    /*WIN32 change by John McMellen*/
	infile=open(filename, O_RDONLY | O_BINARY);

	if (infile==0)
	{
		fprintf(stderr, "*** Error: Cannot open \"%s\"\n", filename);
		return -1;
	}

	n=read(infile,&buffer,2);

	if ((buffer[0]=='P') && (buffer[1]=='K'))
	{
		fprintf(stderr, "*** Error: \"%s\" still appears to be compressed!\n",filename);
		close(infile);
		return -1;
	}

	lseek(infile,0L,SEEK_SET);

//change 07-07-2009 by VE3NCQ to make windows-friendly filenames
	if (ippd==3600)
		sprintf(sdf_filename, "%dx%dx%dx%d-hd.sdf", min_north, max_north, min_west, max_west);
	else
		sprintf(sdf_filename, "%dx%dx%dx%d.sdf", min_north, max_north, min_west, max_west);

//	if (ippd==3600)
//		sprintf(sdf_filename, "%d:%d:%d:%d-hd.sdf", min_north, max_north, min_west, max_west);
//	else
//		sprintf(sdf_filename, "%d:%d:%d:%d.sdf", min_north, max_north, min_west, max_west);

	error=0;
	replacement_flag=0;

	printf("Reading %s... ", filename);
	fflush(stdout);

	for (x=0; (x<=ippd && error==0); x++)
		for (y=0; (y<=ippd && error==0); y++)
		{
			bytes_read=read(infile,&buffer,2);

			if (bytes_read==2)
			{
				if (bil)
				{
					/* "little-endian" structure */

					byte=buffer[0]+(buffer[1]<<8);

					if (buffer[1]&128)
						byte-=0x10000;
				}

				if (hgt)
				{
					/* "big-endian" structure */

					byte=buffer[1]+(buffer[0]<<8);

					if (buffer[0]&128)
						byte-=0x10000;
				}

				/* Flag problem elevations here */

				if (byte<-32768)
					byte=-32768;

				if (byte>32767)
					byte=32767;

				if (byte<=min_elevation)
					replacement_flag=1;

				srtm[x][y]=byte;
			}

			else
				error=1;
		}

	if (error)
	{
		fprintf(stderr,"\n*** Error: Premature EOF detected while reading \"%s\"!  :-(\n",filename);
	}

	close(infile);

	return 0;
}

int LoadSDF_SDF(char *name)
{
	/* This function reads uncompressed
	   SPLAT Data Files (.sdf) into memory. */

	int x, y, n, dummy;
	char sdf_file[255], path_plus_name[255];
	FILE *infile;

	for (x=0; name[x]!='.' && name[x]!=0 && x<250; x++)
		sdf_file[x]=name[x];

	sdf_file[x]='.';
	sdf_file[x+1]='s';
	sdf_file[x+2]='d';
	sdf_file[x+3]='f';
	sdf_file[x+4]=0;

	strncpy(path_plus_name,sdf_path,255);
	strncat(path_plus_name,sdf_file,255);

	infile=fopen(path_plus_name,"rb");

	if (infile==NULL)
		return 0;

	n=fscanf(infile,"%d", &dummy);
	n=fscanf(infile,"%d", &dummy);
	n=fscanf(infile,"%d", &dummy);
	n=fscanf(infile,"%d", &dummy);

	printf("\nReading %s... ",path_plus_name);
	fflush(stdout);

	for (x=0; x<1200; x++)
		for (y=0; y<1200; y++)
			n=fscanf(infile,"%d",&usgs[x][y]);

	fclose(infile);

	return 1;
}

/* BZ2 removed by VE3NCQ to solve compile problems

//char *BZfgets(BZFILE *bzfd, unsigned length)
//{
//	/* This function returns at most one less than 'length' number
//	   of characters from a bz2 compressed file whose file descriptor
//	   is pointed to by *bzfd.  In operation, a buffer is filled with
//	   uncompressed data (size = BZBUFFER), which is then parsed
//	   and doled out as NULL terminated character strings every time
//	   this function is invoked.  A NULL string indicates an EOF
//	   or error condition. */
//
//	static int x, y, nBuf;
//	static char buffer[BZBUFFER+1], output[BZBUFFER+1];
//	char done=0;
//
//	if (opened!=1 && bzerror==BZ_OK)
//	{
//		/* First time through.  Initialize everything! */
//
//		x=0;
//		y=0;
//		nBuf=0;
//		opened=1;
//		output[0]=0;
//	}
//
//	do
//	{
//		if (x==nBuf && bzerror!=BZ_STREAM_END && bzerror==BZ_OK && opened)
//		{
//			/* Uncompress data into a static buffer */
//
//			nBuf=BZ2_bzRead(&bzerror, bzfd, buffer, BZBUFFER);
//			buffer[nBuf]=0;
//			x=0;
//		}
//
//		/* Build a string from buffer contents */
//
//		output[y]=buffer[x];
//
//		if (output[y]=='\n' || output[y]==0 || y==(int)length-1)
//		{
//			output[y+1]=0;
//			done=1;
//			y=0;
//		}
//
//		else
//			y++;
//		x++;
//
//	} while (done==0);
//
//	if (output[0]==0)
//		opened=0;
//
//	return (output);
//}
//
//int LoadSDF_BZ(char *name)
//{
//	/* This function reads .bz2 compressed
//	   SPLAT Data Files into memory. */
//
//	int x, y, dummy;
//	char sdf_file[255], path_plus_name[255];
//	FILE *fd;
//	BZFILE *bzfd;
//
//	for (x=0; name[x]!='.' && name[x]!=0 && x<247; x++)
//		sdf_file[x]=name[x];
//
//	sdf_file[x]='.';
//	sdf_file[x+1]='s';
//	sdf_file[x+2]='d';
//	sdf_file[x+3]='f';
//	sdf_file[x+4]='.';
//	sdf_file[x+5]='b';
//	sdf_file[x+6]='z';
//	sdf_file[x+7]='2';
//	sdf_file[x+8]=0;
//
//	strncpy(path_plus_name,sdf_path,255);
//	strncat(path_plus_name,sdf_file,255);
//
//	fd=fopen(path_plus_name,"rb");
//	bzfd=BZ2_bzReadOpen(&bzerror,fd,0,0,NULL,0);
//
//	if (fd!=NULL && bzerror==BZ_OK)
//	{
//		printf("\nReading %s... ",path_plus_name);
//		fflush(stdout);
//
//		sscanf(BZfgets(bzfd,255),"%d",&dummy);
//		sscanf(BZfgets(bzfd,255),"%d",&dummy);
//		sscanf(BZfgets(bzfd,255),"%d",&dummy);
//		sscanf(BZfgets(bzfd,255),"%d",&dummy);
//	
//		for (x=0; x<1200; x++)
//			for (y=0; y<1200; y++)
//				sscanf(BZfgets(bzfd,20),"%d",&usgs[x][y]);
//
//		fclose(fd);
//
//		BZ2_bzReadClose(&bzerror,bzfd);
//
//		return 1;
//	}
//
//	else
//		return 0;
//}
//
//End VE3NCQ changes */


char LoadSDF(char *name)
{
	/* This function loads the requested SDF file from the filesystem.
	   First, it tries to invoke the LoadSDF_SDF() function to load an
	   uncompressed SDF file (since uncompressed files load slightly
	   faster).  Failing that, it tries to load a compressed SDF file
	   by invoking the LoadSDF_BZ() function. */

	int return_value=-1;

	/* Try to load an uncompressed SDF first. */

	return_value=LoadSDF_SDF(name);

	/* If that fails, try loading a compressed SDF. */

/* Removed by VE3NCQ
	if (return_value==0 || return_value==-1)
		return_value=LoadSDF_BZ(name);
*/

	return return_value;
}

int ReadUSGS()
{
	char usgs_filename[15];

	/* usgs_filename is a minimal filename ("40:41:74:75").
	   Full path and extentions are added later though
	   subsequent function calls. */

//changed 07-07-2009 to make windows-friendly filenames
	sprintf(usgs_filename, "%dx%dx%dx%d", min_north, max_north, min_west, max_west);
//	sprintf(usgs_filename, "%d:%d:%d:%d", min_north, max_north, min_west, max_west);

	return (LoadSDF(usgs_filename));
}

//changed following declaration by VE3NCQ 07-07-2009 to make function legal for my DEV-C++ windows compiler
void average_terrain(int y,int x,int z)
{
	long accum;
	int temp=0, count, bad_value;
	double average;

	bad_value=srtm[y][x];

	accum=0L;
	count=0;

	if (y>=2)
	{
		temp=srtm[y-1][x];

		if (temp>bad_value)
		{
			accum+=temp;
			count++;
		}
	}

	if (y<=mpi)
	{
		temp=srtm[y+1][x];

		if (temp>bad_value)
		{
			accum+=temp;
			count++;
		}
	}

	if ((y>=2) && (x<=(mpi-1)))
	{
		temp=srtm[y-1][x+1];

		if (temp>bad_value)
		{
			accum+=temp;
			count++;
		}
	}

	if (x<=(mpi-1))
	{
		temp=srtm[y][x+1];

		if (temp>bad_value)
		{
			accum+=temp;
			count++;
		}
	}

	if ((x<=(mpi-1)) && (y<=mpi))
	{
		temp=srtm[y+1][x+1];

		if (temp>bad_value)
		{
			accum+=temp;
			count++;
		}
	}

	if ((x>=1) && (y>=2)) 
	{
		temp=srtm[y-1][x-1];

		if (temp>bad_value)
		{
			accum+=temp;
			count++;
		}
	}

	if (x>=1)
	{
		temp=srtm[y][x-1];

		if (temp>bad_value)
		{
			accum+=temp;
			count++;
		}
	}

	if ((y<=mpi) && (x>=1))
	{
		temp=srtm[y+1][x-1];

		if (temp>bad_value)
		{
			accum+=temp;
			count++;
		}
	}

	if (count!=0)
	{
		average=(((double)accum)/((double)count));
		temp=(int)(average+0.5);
	}

	if (temp>min_elevation)
		srtm[y][x]=temp;
	else
		srtm[y][x]=min_elevation;
}

void WriteSDF(char *filename)
{
	/* Like the HGT files, the extreme southwest corner
	 * provides the point of reference for the SDF file.
	 * The SDF file extends from min_north degrees to
	 * the south to min_north+(mpi/ippd) degrees to
	 * the north, and max_west degrees to the west to
	 * max_west-(mpi/ippd) degrees to the east.  The
	 * overlapping edge redundancy present in the HGT
	 * and earlier USGS files is not necessary, nor
	 * is it present in SDF files.
	 */

	int x, y, byte, last_good_byte=0;
	FILE *outfile;

	printf("\nWriting %s... ", filename);
	fflush(stdout);

	outfile=fopen(filename,"wb");

	fprintf(outfile, "%d\n%d\n%d\n%d\n", max_west, min_north, min_west, max_north);

	for (y=ippd; y>=1; y--)		/* Omit the northern most edge */
		for (x=mpi; x>=0; x--) /* Omit the eastern most edge  */
		{
			byte=srtm[y][x];

			if (byte>min_elevation)
				last_good_byte=byte;

			if (byte<min_elevation)
			{
				if (merge)
				{
					if (ippd==3600)
						fprintf(outfile,"%d\n",usgs[1200-(y/3)][1199-(x/3)]);
					else
						fprintf(outfile,"%d\n",usgs[1200-y][1199-x]);
				}

				else
				{
					average_terrain(y,x,last_good_byte);
					fprintf(outfile,"%d\n",srtm[y][x]);
				}
			}
			else
				fprintf(outfile,"%d\n",byte);
		}

	printf("Done!\n");

	fclose(outfile);
}

int main(int argc, char **argv)
{
	int x, y, z=0;
	char *env=NULL, string[255], *s=NULL;
	FILE *fd;

	if (strstr(argv[0], "srtm2sdf-hd")!=NULL)
	{
		ippd=3600;	/* High Definition (1 arc-sec) Mode */
		strncpy(string,"srtm2sdf-hd\0",12);
	}

	else
	{
		ippd=1200;	/* Standard Definition (3 arc-sec) Mode */
		strncpy(string,"srtm2sdf\0",9);
	}

	mpi=ippd-1;		/* Maximum pixel index per degree */

	if (argc==1 || (argc==2 && strncmp(argv[1],"-h",2)==0))
	{
		if (ippd==1200)
			fprintf(stderr, "\nsrtm2sdf: Generates SPLAT! elevation data files from unzipped\nSRTM-3 elevation data files, and replaces SRTM data voids with\nelevation data from older usgs2sdf derived SDF files.\n\n");

		if (ippd==3600)
			fprintf(stderr, "\nsrtm2sdf-hd: Generates SPLAT! elevation data files from unzipped\nSRTM-1 elevation data files, and replaces SRTM data voids with\naverages, or elevation data from older usgs2sdf derived SDF files.\n\n");

		fprintf(stderr, "\tAvailable Options...\n\n");
		fprintf(stderr, "\t-d directory path of usgs2sdf derived SDF files\n\t    (overrides path in ~/.splat_path file)\n\n");
		fprintf(stderr, "\t-n elevation limit (in meters) below which SRTM data is\n\t    replaced by USGS-derived .sdf data (default = 0 meters)\n\n");
		fprintf(stderr, "Examples: %s N40W074.hgt\n",string);
		fprintf(stderr, "          %s -d /cdrom/sdf N40W074.hgt\n",string);
		fprintf(stderr, "          %s -d /dev/null N40W074.hgt  (prevents data replacement)\n",string);
		fprintf(stderr, "          %s -n -5 N40W074.hgt\n\n",string);

		return 1;
	}

	y=argc-1;

	min_elevation=0;

	for (x=1, z=0; x<=y; x++)
	{
		if (strcmp(argv[x],"-d")==0)
		{
			z=x+1;

			if (z<=y && argv[z][0] && argv[z][0]!='-')
				strncpy(sdf_path,argv[z],253);
		}

		if (strcmp(argv[x],"-n")==0)
		{
			z=x+1;

			if (z<=y && argv[z][0])
			{
				sscanf(argv[z],"%d",&min_elevation);

				if (min_elevation<-32767)
					min_elevation=0;
			}			 
		}
	}

	/* If no SDF path was specified on the command line (-d), check
	   for a path specified in the $HOME/.splat_path file.  If the
	   file is not found, then sdf_path[] remains NULL, and a data
	   merge will not be attempted if voids are found in the SRTM file. */

	if (sdf_path[0]==0)
	{
//changed by VE3NCQ 07-07-2009 to make windows friendly
		char* hpath;
        env=getenv("HOMEPATH");
		hpath=getenv("HOMEDRIVE");
		sprintf(string,"%s%s\\.splat_path",hpath,env);
//		env=getenv("HOME");
//
//		sprintf(string,"%s/.splat_path",env);

		fd=fopen(string,"r");

		if (fd!=NULL)
		{
			s=fgets(string,253,fd);

			/* Remove <CR> and/or <LF> from string */

			for (x=0; string[x]!=13 && string[x]!=10 && string[x]!=0 && x<253; x++);
			string[x]=0;

			strncpy(sdf_path,string,253);

			fclose(fd);
		}
	}

	/* Ensure a trailing '/' is present in sdf_path */
//change by VE3NCQ 07-07-2009 windows uses backslash \ (must be escaped)
	if (sdf_path[0])
	{
		x=strlen(sdf_path);

		if (sdf_path[x-1]!='\\' && x!=0)
		{
			sdf_path[x]='\\';
			sdf_path[x+1]=0;
		}
	}

	if (ReadSRTM(argv[z+1])==0)
	{
		if (replacement_flag && sdf_path[0])
			merge=ReadUSGS();

		WriteSDF(sdf_filename);
	}

	return 0;
}

